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SUMMARY OF RESEARCH RESULTS 

1. The following report titled "Computation of Optimal Low- and Medium- Thrust 
Orbit Transfers" gives the detail of the research results. We summary these 
results and future research plan here. 

2. The first-order necessary conditions for a general final mass maximization 
problem has been set up. In the problem formulation we include second- 
harmonic oblateness, atmospheric drag, and allow three-dimensional, non- 
coplanar, non-aligned elliptic orbits. In order to ease the numerical calculation 
we transform the original free final-time problem to a fixed final-time problem, 
and non-dimensionalize the state variables. 

3. Although we can use the constant angular momentum equation, the 
conservative energy equation, and the orbit equation to specify the boundary 
conditions for the terminal orbit, we notice that this set of boundary conditions 
does not uniquely determine an orbit. This is due to the fact that for a given 
point in space we can have two different velocity vectors (difference in direction 
only) and yet have the same angular momentum and energy. Proper boundary 
conditions should be three eccentricity vector equations plus three angular 
momentum vector equations. Since both eccentricity and angular momentum 
equations specify the same orbit plane, one of these equations is redundant. That 
is for a three dimensional problem we only need five equations out of both sets 
of equations. For two dimensional problem we need two eccentricity equations 
and one angular momentum equation. 

4. We have applied two indirect optimization methods: BOUNDSCO and MBCM 
(minimizing-boundary-condition method) successfully to several simplified 
examples. The examples are two dimensional with oblateness effect and 
atmospheric drag force. Both methods converge to the solutions with about the 
same sensitivity in the initial guess. Although we have more freedom in 
selecting the initial guess at every node points, BOUNDSCO does not adjust the 
number of switching points and the switching pattern during the iteration. On 
the other hand, MBCM implements the switching function into the integrator 
and adjust the switching points and the switching pattern automatically during 

the iteration. 

5. Our current plan is to combine advantageous features of BOUNDSCO and MBCM 
into a new algorithm. The new algorithm will use the idea of the multiple-point 
shooting method to spread the unknowns among the node points, and between 
two node points applies the minimizing-boundary-condition method. 

6. There is still a question about the local optimum or global optimum for free final 
time problem. We have some difficulty in converging the transversality 
condition for the free final time case. In Edelbaum's paper, he shows that three 
impulses control is usually minimum. However, such claim for low and 
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medium thrust has not been shown anywhere. Our current hypothesis suggests 
S lobal solution will be at infinite final time and local 

minimum solutions exist for finite final time. We expect to answer this 
question by obtaining all the local minimum solutions (if they exist) and 
compare their cost functions along the final time axis. 
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Computation of Optimal Low- and Medium- 
Thrust Orbit Transfers 


Abstract 

This report presents the formulation of the optimal low- and medium- 
thrust orbit transfer control problem and methods for numerical solution of 
the problem. The problem formulation is for final mass maximization and 
allows for second-harmonic oblateness, atmospheric drag, and three- 
dimensional, non-coplanar, non-aligned elliptic terminal orbits. We setup 
some examples to demonstrate the ability of two indirect methods to solve 
the resulting TPBVPs. 

The methods demonstrated are the multiple-point shooting method as 
formulated in H. J. Oberle's subroutine BOUNDSCO, and the minimizing 
boundary-condition method (MBCM). We find that although both methods 
can converge solutions, there are trade-offs to using either method. 
BOUNDSCO has very poor convergence for guesses that do no exhibit the 
correct switching structure. MBCM, however, converges for a wider range of 
guesses. However, BOUNDSCO's multi-point structure allows more freedom 
in guesses by increasing the node points as opposed to only guessing the 
initial state in MBCM. Finally, we note an additional drawback for 
BOUNDSCO: the routine does not supply information to the users routines 
for switching function polarity but only the location of a preset number of 
switching points. 


I. INTRODUCTION 


The ability to perform any given orbit transfer with a minimum use of 
fuel is obviously desirable. Useful solutions to this problem will account for 
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at least some approximation to real-life. Therefore, a formulation that 
indudes second-harmonic oblateness and atmospheric drag will be useful. 

This report follows such a derivation all the way through to the 
establishment of a two-point boundary-value problem for optimal low- and 
medium- thrust orbit transfer. The core cost function is defined simply as the 
final mass of the spacecraft plus fuel, setting the tone for the maximization 
problem. The differential constraint is thoroughly defined in terms of the 
oblateness model and an assumed atmosphere model. 

The thrust (control) appears linear in the differential constraint. This 
results in bang-bang control or singular-arc solutions for the final mass 
maximization problem. Although bang-bang control is assumed here the 
possibility of having a singular arc has not been ruled out for a general case. 
In order to ensure the singular arc solution does not occur, we check the 
derivative of the switching function at each switching point. However, when 
our programs reach a non-optimal solution high frequency chattering 
solutions do occur occasionally. This could indicate that singular-arc 
solutions are possible for some modification of system parameters and 
models. 

The final mass maximization problem should be a free final time 
optimal control problem. For impulsive thrust, the Hohmann transfer gives 
minimum fuel but maximum transfer time. Although the three-impulse 
Hohmann transfer performs better than a two-impulse Hohmann transfer, 
Edelbaum 1 shows that the number of impulses may be finite for a global 
minimum, for low- and medium-thrust orbit transfer the same conclusion 
has not be shown anywhere. One hypothesis is that the global minimum will 
be at infinite final time and local minimum solutions exist for finite final 
time. In other words, this assumes for a given number of switching points 
(must be at least two) there is a local minimum with finite final time. We do 
have difficulties in converging the transversality condition corresponding to 
optimal final time. 

We present solutions to three specific optimization problems. These 
solutions represent the ability of the two TPBVP solvers. The methods 
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considered are (1) BOUNDSCO, a multi-point shooting algorithm devised by 
H. J. Oberle and (2) the minimizing boundary-condition method (MBCM), a 
modification to the shooting method devised by the authors of ref(9). 

Both methods converge solutions for about equal sensitivity in initial 
guesses. In order to achieve the same accuracy along the path, BOUNDSCO 
needs to converge the boundary conditions at every node point to the same 
accuracy as the integration routine, the number of switching points and the 
switching pattern need to be assumed and stay unchanged when BOUNDSCO 
is used. On the other hand, MBCM does not constrain the number of 
switching points and MBCM updates the switching pattern along the 
integrated path. 


II. THE PROBLEM 

The problem discussed herein is the following: maximize the final 
mass of a thrusting spacecraft for a given orbital transfer. The craft can be 
considered as under the influence of some planet's gravitational field and 
atmospheric drag. The thrust of the spacecraft is limited between zero and 
some Tmax- The transfer will be defined by the terminal orbits. Solutions are 
sought for both fixed and free final time problems and both the case of fixed 
and free terminal points. 

n. 1 . The Cost Functional 

The core cost functional must be defined. We shall define the cost as 

J = m(tf) ' (1) 

where m(tf) represents the' mass of the spacecraft plus its fuel at the end of the 
orbital transfer. We shall use the methods of optimal control to maximize 
the cost functional, thereby maximizing the final mass and solving the 
problem. 
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n. 2. Differential constraint: System Dynamics 

We represent the spacecraft by a point mass and assume it to be a 
thrusting craft acted upon the by the aerodynamic drag and oblate-body 
gravity forces of a central body. We also represent the central body, or planet, 
as a point mass positioned at its own Center of gravity. We restrict the 
problem to crafts of mass much smaller than that of the central body, 
allowing us to fix the planet in inertial space. We shall describe this inertial 
space with a rectangular Cartesian inertial reference frame. All motion 
within this frame of reference agreeing with the above assumptions must 
satisfy the well-known Newton's equation: 


jj_d(rnv) (2) 

dt 

where m is the spacecraft mass and v is its velocity with respect to the 
reference frame. 


In this case, gravity, drag, and thrust make up the total force acting on 
the craft. The thrust on the craft is composed of two separated thrusts, the 
pressure thrust and the thrust created by the expulsion of mass. That is. 


F = Fpressure thrust " F*ag “ Fg rav ity — mv - lfav e 
where is the expulsive velocity of mass. Therefore, 

F total thrust s Fpressure thrust + mv e 

and 

mV = F thjust ■ F<Jrag * Fg ra vity 


( 3 ) 


( 4 ) 


( 5 ) 


We write the total thrust, herein referred to as just thrust, as some 
time-varying magnitude, T, independent of a time-varying direction, e: 



5 


F thrust — T C (6) 

Note that e is expressed as a unit vector. For a three dimensional thrust 
vector the control requires three components. For two dimensional problems 
only two independent control components are required. 

The mass will decrease according to 


m 


. T. 

go lsp 


(7) 


We assume that the atmosphere surrounding the central body can be 
described by an exponential model of the standard atmosphere 2 . The 
following equation 3 describes such a drag force: 


Fdrag^Poe-^^SCDfl 2 ! = lp.c-^-^SCbHv w 

where fl and r Q are constants from the atmosphere model describing air 
density variation in the prescribed altitude region, p G is the atmosphere 
density for the altitude r 0 , S is the wetted area of the craft, Cd is the craft s drag 
coefficient, and v is craft's current velocity with respect to the inertial 
reference frame. We are assuming that no matter the orientation of the craft 
the product of SCd remains the same and that the craft always remains in a 
region where the chosen exponential atmosphere model is valid. 

Within the confines of this study, the only other influence on the craft 
is gravitational potential energy. The gravitational potential energy to the 
second harmonic is 4 : 

u= .5^..1jr^(1-3co S 29) (9) 

Where R is the equatorial radius of the central body, 6 is the latitude angle of 
the curreo position from the equator, and r is the distance from the cential 
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body's center of gravity to the current position of the craft with respect to the 
inertial reference frame, p is the gravitational constant for the central body, m 
is the mass of spacecraft, and J is a constant describing the mass distribution of 
the planet. There are additional mass distribution terms but we shall truncate 
the series here. 

I ; 

We now assume that (1) the central body is fixed at the center of the 
reference frame and (2) that the plane of the central body's equator is aligned 
with the x-y plane. The assumption (1) means that the position, velocity, and 
acceleration of the craft are now measured with respect to the central body. 
The assumption (2) means that we may describe r with Cartesian coordinates 
by 


r = V x 2 + y 2 + z 2 

and we may describe 0 with Cartesian coordinates by 

z = r cos 0 

We may now write the gravitation potential as 

^ r 


( 10 ) 


( 11 ) 


( 12 ) 


au 


The force experienced by moving in a potential field U is yy 
Performing this operation on the gravitational potential yields 


gravity 


au\T 
a? 


( 13 ) 


and 
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( 14 ) 


3U _m(i 

dy 


^ ytJ R^yO - 5 ©■) 


(15) 




(16) 


All of the dynamics combines to form the following equations of motion 


m x = Te, • x - J R 2 ^x (1-5 Iff ) - ^ p» e P(r ,b) SCdvx 


my = Te y -^y - JR 2 ^y (1 - 5(f) 2 ) - Ip^e^SCovy 
mz = Te, -^z - JR 2 ^z( 3 - 5(^) - Ip^e ^SCdvz 

«3 *5 'i' ? 


(17a) 

(17b) 

(17c) 


which can be written in vector-matrix form as 




f 

1 

1 

t-rt 

0 

0 


t ■* 
r = — e 

- \Lt - , 

JR 2 ^ 

0 

>- s (f] 

« 

l 

m 

r 5 

r 5 

J 




0 

0 

3 -- 5 ft 

rJI 


i^e'^-'c’SCovr (18) 
2 m 


To conform to convention we make the change from J to J 2 as described in 
ref(4): 



r 


( 19 ) 




1 0 o' 
0 1 0 
VL 0 0 3 - 



lfi. e -P(r -*.>sC DV r 
2 m 


This can also be rewritten as a first-order system: 


(20a) 






1 0 O' 
0 1 0 
VL o o 3. 



- i^e-^SCovv 
2 m 


(20b) 


m. THE FIRST-ORDER necessary conditions 

All problems herein perform a maximization of the final mass. Now, 
in order to write the adjoined cost functional we need to know what is 
included in the state, in the control, and what the constraints on these are. 

First, however, we note that the problems herein are also free-final- 
time problems. The three differential equations above are written with 
respect to the independent variable t (time). For ease in numerical methods, 
we want to transfer the problem from free- to fixed- final time. This means 
that we must define a new independent variable x (non-dimensional time) to 
be used in the place of t (dimensional time). This allows tf to become a state 
variable. We make the following scaling: 

t = t f t . (2D 

Therefore, to translate this to a fixed-time problem, tfmust multiply the 
derivatives of the states. The dot above a variable now means a derivative 
with respect to t. 


We know what to include in the state, x(x): 
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x(T)=[ F(t) V(x) m W tf 7 (22) 

We also know that our state is confined by the system dynamics so that 


X(T) =tf 



f 1 0 0 ' 
0 1 0 

' 5 

\ 2 JS \ 

l 0 0 3 - 

KrIJf 


5 T - e- S C D 1 


.JL- 


£oIq> 

0 


(23) 


for all time x e [0,1]. This is the differential constraint of the control problem. 

The thrust magnitude has both an upper and a lower bound. The 
upper bound we shall call T m ax/ the lower bound is obviously zero. We, 
therefore, also have an inequality constraint that must be satisfied for all time 
x e [0,1]: 


(T - T max )T < 0 

and Eqn (25) can be rewritten as an equality constraint 


(24) 


(T - T max )T + a 2 = 0 (25) 

where a is a slack variable, free to change with time. Finally, we need to 
specify the terminal orbits. We will do so by writing a vector equation 


y(x(0),x(l)) = 0 


(26) 
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that is only satisfied when our initial and final states both lie on their 
respective orbits. 

Now we know enough to write the adjoined cost functional for this 
problem: 


J = m(l) + \mj/(x(0),x(1)) + 



t),u(t)) - x] + m[(T - T max )T 



dx 


(27) 


where m(l) is the final mass and y(x(0),x(l)) = 0 represents the boundary 
conditions. 

The X shown in the cost functional is the costate vector, also called the 
Lagrange multipliers. This vector will be of the same dimensions as the state. 
For simplification's sake, we will segment this vector as follows: 


X(t) = 


— T 

K (t) 


Xv (t) ^m(0 



(28) 


Also, in Eqn (27), v is the Lagrange multiplier corresponding to the boundary 
conditions vy. v is a constant in time. 


m.l. The Hamiltonian 


With the pertinent dynamics defined, we are now able to write the 
Hamiltonian for the system. We take the Hamiltonian from the cost 
functional as 

_rr 

H(x(t),u(t)) = X [f(x(t),u(t))] + p[(T - T m ax)T + a 2 ] (29) 


A major simplification can be made now. Notice that, excluding the 
constraint on the thrust, the Hamiltonian is linear with respect to the control 
T (but we shall see it is not linear with respect to e ): 
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H(x(t),u(t)) = e - 


(30) 


This, in conjunction with the structure of the thrust constraint, means 
that we may assume that this is a 'bang-bang' control problem. Enough is 
known about this type of problem so that we may do the following: 

(1) Define a new Hamiltonian that differs from the original only in the 
omission of the thrust constraint. 

H(x(t),u(t)) = X [f(x(t),u(t)>] ^ 


(2) Establish what will be called the switching function. In general, the 
switching function is defined by the partial derivative of the Hamiltonian 
with respect to the control by which it is linear. For this problem. This is 

3H 

done by evaluating 




Isp go 


(32) 


This, Eqn (32), is the switching function. 

(3) Evaluate a restricted case of the well-known Euler-Lagrange 
equations. Most of these determine the costate dynamics and we shall see 
these in section III.2, however, the last one determines part of the control for 
the problem. This equation is 



(33) 


Evaluating (33), it appears that H is linear with respect to e. However, we 
must remember that e represents only the direction of the thrust. If we 
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exchange e for some angle 6 and define this angle as between e and X v we 
may write 


H 

H 


m 


Xv 


-I-lxJcosQ + 
m r 


l k 

cos 0 + • • • 

(34) 

since 

H - 1 

(35) 


Evaluating He we find that 


— = -^-0si n0 (36) 

30 m 

and this equals zero only when the vectors are parallel. There are only two 

choices for e: in the direction of K or in the exact opposite direction. Since we 
are maximizing the final vehicle mass, we need to have H e0 negative (one of 
the sufficient conditions for the second variation). This is only satisfied with 

e in the direction of X v or 


(37) 

We must obey this for all time t e [0,1]. This result is consistent with that of 
Lawden's primer vector 5 . 

(4) Perform bang-bang control with T. This means that T is always on- 
boundary, i.e. T=0 or T=T max at any t e 10,11- We know which value to use 
for T by evaluating the switching function, which we can now write as 

„ Xn, (38) 

Ht - m- ' ^ 



The bang-bang control law is 



4 


(39) 
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H t £ 0 


T = X 


max 


Hj < 0 T = 0 


This switching structure satisfies the Pontryagin maximum principle by 
maximizing the Hamiltonian using T. .■ ; 


m.2. The Costates 

The costate dynamics can be found from the following Euler-Lagrange 
equations, relating them to the Hamiltonian: 


K = - 


dH T 

a? 


(40) 


X v 


3HT 

9v 


(41) 


K = - 


8H 

dm 


(42) 


To evaluate these, we must first substitute the equations of motion into the 
Hamiltonian 


H = tflTv 




( 43 ) 


+ 


_x_ 

go lsp 


) 
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When evaluated, these become the following vector and scalar differential 
equations: 






- ***>SCdv 

2 m r , 




- [ °' 
where k = 0 

. 1 . 


Xy 


= % 



1 £i g-P (r-r 0 ) 

2 m 



( 45 ) 


X 


m 


= tf 



— T _ 

Xy 0 + 


liLe-PNisCbv xj\ 

2 m 2 


( 46 ) 


IV. SOLVED PROBLEMS 

IV.l. Simplifications 

We have made a few simplifications that ease the formulation of the 
numerical problem and its solution. The' first of these is the reduction from a 
three-dimensional problem to a two-dimensional problem. To remove this 
dimension, we simply remove the z-component to all equations. Because of 
the chosen coordinate system, this also means that all orbit transfers 
considered are equatorial. Unfortunately, the effect of oblateness is 
substantially decreased for this case. The other simplification is the restriction 
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of problems to fixed-initial points. This also greatly eases the problem 
formulation. The third and final simplification is the fixing of the final time. 

IV.2. The Two-Point Boundary Value Problem 

As a result of the simplifications, the boundary conditions have been 
stated in two dimensions. The starting orbit determines the initial conditions 
on position and velocity. The final conditions, however, require a more 
abstract specification as we do not know exactly at what point the craft will 
enter this orbit. The following relations specify the final orbit: (All of the 
following conditions is to be evaluated at the final time, tf, or t=l.) 

R x V = H | 

(Angular Momentum) yl: (<x, y> x <u, v> = xv - yu = h: (47a) 

xv - yu - h = 0 | 


e x = -M(v 2 -?)x - (f^vju] 

r 

(eccentricity vector (x) ) y2: ( (47b) 

b-[( v2 -?) x • - e * = 0 


•» - ^[(v 2 -f)y - (^H 

(eccentricity vector (y) ) y3: ( ' j (47c) 

ut( v2 '?) y ‘ " “y = 0 

I r 


Note that the orbit equation for x-axis aligned orbits and the energy 
equation can replace (47b) and (47c). However, the combined constraints of 
angular, momentum, orbit, and energy equations do not uniquely specify an 
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' x-axis aligned final orbit. There are two possible velocity vectors at one point 
with the same angular momentum and energy. 

These conditions completely determine the final orbit. However, these 
conditions do not complete the two-point boundary-value problem. To 
complete the TPBVP, the methods of optimal control supply use with a set of 
natural boundary conditions found by evaluating 


I(1) ■ (sfj 


(48) 


where G is the constructed from the function portion of the cost functional, 
e.g. for the cost functional 


J = m(l) + \Mj/(x(0),x(l)) + 



(49) 


Gis 


G = m(l) + \mi/(x(0),x(1)) 


(50) 


Constructing G with the above conditions on the states, we can find 
conditions on the costates at t=l: 


G = m + Vi (x v - y» - h) +[ v 2 v 3 ]Mv 2 -f)f - (?v)v) - e‘ 
evaluating Eqn(49) gives 


(51) 




3G 

3x 


Vi (v) + v 2 





& 


+ 


U V 


(52a) 



9G 

3y 


= Vi (-u) + v 2 




(52b) 
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*U = = Vi(-y) + v 2 j 1 |pj + V 3 ^(2yu-xv) 

= ^ = Vi(x) + v 2 ^(2 x v * u y) + v 3 |^j 


(52c) 


(52d) 


(52e) 


note that the constant Lagrange multipliers Vi are additional unknown's. 

The last condition deals with the final time. If the final time were free 
we would use the transversality condition 


H(x(l),u(lU(l)) = 


3G 

dtf 


(52f) 


or, for this problem 


H(x(l),u(l),Ml))= 0 

However, all the solutions presented in this report are fixed-final time. Note, 
however, that the same algorithm can be used for both types of problems, all 
that is required is that equation (53g) be replaced by the specification of tf. 


IV.3. Non-Dimensionalization ' 

To improve accuracy, we have non-dimensionalized the problem. 
This aids in a few ways. First, the integration of the state is more accurate 
because all variations are on the same order. Second, convergence is 
improved because all the boundary conditions are immediately placed at or 
near the same order. Our non-dimensionalized parameters are as follows: 
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r 



(53a) 



ms-ffi- 

m* 


(53b) 


(53c) 


fr _ tf (53d) 

W 

and they require the following 


T = T / 

(53e) 

[ i/r* 2 


(go Isp) - (go Isp) / \J^ 

(53f) 

r = 1° 

r ° ~r* 

(53g) 

P = P r * 

(53h) 

((V.CdSJeIp.CdS)^- 

(53i) 

R=^ 

r w 

(53 j) 
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The choices of r* and m* are completely arbitrary. However, it needs to be said 
that after a problem is solved by these nondimensionalizations, rescaling 
must be excersized with caution. This is a direct result of the atmosphere 
model; if the rescaling is not consistent with the atmosphere model, the 
results are invalid, e.g. rescaling also rescales the atmosphere model (note Eqn 
( 8 )). 


If we solve Eqs (53a-j) such that the dimensional parameter is on the 
left-hand side and then substitute into the original dynamics we find 
equations that are exactly equal to the original equations with p=l (The value 
of J 2 , however, has no dimensions and is not changed). This can be extended 
to the boundary equations and the costate differential equations. A special 
note is required for the costates: the costates resulting from the solution to the 
problem with this transformation will be some scalar multiplied by the 
'dimensional' costates, e.g. 



(53k) 


which requires 


v 



(531) 


where X* is completely arbitrary. This is easily verified by substitution into 
the differential equations and boundary conditions. 


IV.4. Atmosphere Model 

Any atmosphere is usable by simple substitution early in the 
derivation of the differential constraint. For the purposes of this report we 
have chosen a very simple atmosphere model. The model is not intended to 
accurately represent the Earth's atmosphere, or any other planet for that 
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matter. It is implemented only for the purpose of demonstrating the 
methods for solving the optimization problem. 

Our model is defined at 450km altitude above the planet's equator. The 
entire atmosphere is assumed isothermal. The temperature is 1000K. The 
density at the definition altitude is 1.184xl0' 12 kg/m 3 . This definition point 
for this model is taken from the 1976 U.S. Standard Atmosphere 6 . The 
atmosphere is assumed spherical above the oblate planet. For real-world 
solutions, we strongly recommend the use of the latest standard atmosphere 
or some appropriate approximation thereof. The contemporary standard 
atmosphere can be found in ref (7). 

IV.5. The Multiple Point Shooting Method of BOUNDSCO 

One method we are currently using to attempt to solve the TPBVP is 
the multiple point shooting method. The specific algorithms we are 
currently using are those given by H. J. Oberle in his subroutine 
BOUNDSCO 7 , written in FORTRAN. His method, a complete description of 
which can be found in ref(8), is a modification of the traditional well-known 
multiple point shooting method. 

The use of this method requires the writing of a few routines that 
define the problem. These routines include, of course, the calling program 
itself, a subroutine defining the differential constraint (or system dynamics), 
and a subroutine that defines the constraints on the problem. 

The state used in BOUNDSCO differ slightly from the state defined in 
this report. We have simply adjoined the v vector to the state. This requires 
also that the system dynamics included a corresponding number of zero 
derivatives. We justify for this by noting that it allows the statement of the 
absolute and natural boundary conditions exactly as they are in this report. If 
we did not implement this, we would have to solve the system of three of the 
natural boundary conditions for the v and substitute the result into the 
fourth equation, using it in place of the four. This may seem desirable, one 
equation in the place of four, however, the simple structure of the four 
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equations is much more desirable than the complex structure of the one 
equation. 

There is one particular feature that makes BOUNDSCO attractive: the 
explicit inclusion of switching points in the problem formulation. Oberle 
allows the user to specify the switching, function outside of the system 
dynamics. This simplifies integration and improves convergence. There is a 
tradeoff; the user must assume a switching structure and verify it outside of 
BOUNDSCO. 

IV.6. The Minimizing-Boundary-Condition Method 

The second method we are using is called the Minimizing-Boundary- 
Condition Method (MBCM) 8 . It is described in ref(9). This method is a 
modification to the shooting method. It expands the set of available solutions 
by removing one boundary-condition. The choice of this boundary-condition 
is arbitrary. Since there is a much larger set of solutions, it is much easier to 
solve the resulting boundary-value problem. Once this is accomplished, the 
search for the solution that incorporates the final boundary conditions is 
treated as a minimization problem. The gradient is numerically calculated 
and used to update the initial state until the last boundary condition is 
satisfied. This method is at least as effective as BOUNDSCO in solving the 
two-point boundary-value problems for the current solved optimal orbit 
transfers. 

The switching structure of optimal control is included in MBCM. The 
program checks the switching function at each integration step. If the 
switching function alters sign at one integration step, the program stops the 
integration and restores all the states to 'the beginning of the step. A secant 
method then calculates a smaller step size for integrating the switching 
function to an exact zero point. From our experience with MBCM some 
sensitive problems need fourteen digits of accuracy in their switching 
function. Once the integration passes the switching points the program 
switches the control and uses a normal step size for integration. 


IV.7. Sample Problems and Solutions 
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Several solutions are presented in this section, all of which both 
methods were able to converge. As a matter of fact, in most cases, the 
solution to one problem can be used as the guess to a different problem and 
the program(s) will converge. All problems have been nondimensionalized 
and use the atmosphere model presented above. 

The first problem presented is a fixed-final-time cirde-to-circle orbit 
transfer: 

Find an extremal for the maximum final-mass problem which travels from a 
circular orbit of a=3.847 at y=3. 72 to another circular orbit of a-1.5. 
The available thrust is (a) 0.9, (b) 02 and goIsp=51 254. The initial 
mass is 1527. The allowed time for transfer is 125. p 0 SCo=3. 894x10' 

17 

The optimal trajectories are shown in Fig. 1 for T=0.9 and in Fig. 3 for 
T=0.2. Their switching functions are shown in Fig. 2 for T=0.9 and in Fig. 4 
for T=0.2. 


The second problem presented is a fixed-final-time apse-aligned ellipse- 
to-ellipse orbit transfer: 

Find an extremal for the maximum final-mass problem which travels from an 
orbit of a=3.847 and rp— 3.756 at y=3.76 to another orbit of a=15 and 
Tp-1. The apses of the orbits are aligned with the x-axis. The 
available thrust is (a) 0.9, (b) 02 and g 0 I S p=51.254. The initial mass 
is 1527. The allowed time for transfer is 12. PoSCd- 3. 894x1 0' J7 . 

The optimal trajectories are shown in Fig. 5 for T=0.9 and in Fig. 7 for 
T=0.2. Their switching functions are shown in Fig. 6 for T=0.9 and in Fig. 8 
for T=0.2. 

The third problem presented is a fixed-final-time non-apse-aligned 
ellipse- to-ellipse orbit transfer: 
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Find an extremal for the maximum final-mass problem which travels from an 
orbit of a=3.847 and r p =3.756 at y=3. 76 to another orbit of a=l _5 and 
r p =l. The apses of the initial orbit is at an angle of 153° with the x- 
axis, clockwise. The apse of the final orbit is at an angle of 109 0 with 
the x-axis, counter-clockwise. The available thrust is (a) 0.9, (b) 02 
and g 0 l sp-5 1254. The initial' mass is 1J27. The allowed time for 
transfer is 10. PoSCd=3.894x10- 17 . 

The optimal trajectories are shown in Fig. 9 for T=0.9 and in Fig. 11 for 
T=0.2. Their switching functions are shown in Fig. 10 for T=0.9 and in Fig. 12 
for T=0.2. 

V. CONCLUSIONS 


The performance of BOUNDSCO was mixed. The ability of the routine 
to converge solutions is quite strong, however, there is a flaw. BOUNDSCO 
does not supply information to the user's routine concerning the polarity of 
the switching function. The user must assume in all his/her code that the 
desired switching structure is correct. The result of this is that BOUNDSCO 
often allows itself to converge solutions with inconsistent switching 
functions. This would not be so bad, except for one other difficulty with 
BOUNDSCO: the routine does not attempt to aid the user in any way with the 
initial guess. For example, one finds it nearly impossible to converge a two- 
burn solution without the insight to guess an initial state that, when 
integrated, produces two crossings of the switching function (this is actually, 
not too difficult, if one pays attention to the sign of the switching function 
and its derivative when making guesses). However, when BOUNDSCO does 
produce correct solutions, they are as accurate as the user can specify. The 
solutions presented above satisfy their boundary conditions within 10‘ 14 
absolute error. 


The performance of the minimizing-boundary condition method was 
also quite promising. This method has one distinct advantage over 
BOUNDSCO, it explicitly disallows inconsistent switching functions. The 
method checks the switching function during, but separately from, 
integration to determine where the switching points are and, most 
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importantly, what the switching function polarity is. This method is, 
however, currently a simple shooting method and it exhibits the difficulty of 
the same. It is expected that if the method is extended to a multiple-point 
shooting method, it's performance will rival, if not exceed that of 
BOUNDSCO. 

»« •. 

And thereby we come to the recommendation of this study: the 
development of a method that is a hybrid of multiple-shooting and the 
minimizing-boundary-condition method. 
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Fig. 5 Mass-Optimal Aligned-Ellipse-to-Ellipse Orbit Transfer, T=0.9 



Fig. 6 Graph of Switching Function Corresponding to Figure 5, T=0.9 
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Fig. 7 Mass-Optimal Aligned-Ellipse-to-Ellipse Orbit Transfer, T-0.2 



Fig. 3 Graph of Switching Function Corresponding to Figure 7, T=0.2 
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Fig. 9 Mass-Optimal Non-Aligned-Ellipse-to-Ellipse Orbit Transfer, T=0.9 
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Fig. 10 Graph of Switching Function Corresponding to Figure 9, T=0.9 









X coordinate 


Fig. 11 Mass-Optimal Non-Aligned-Ellipse-to-Ellipse Orbit Transfer, T=0.2 









